USE ElliDef
USE NodeInfoDef
USE metric

USE Elliptic_Solvers
IMPLICIT NONE
TYPE (NodeInfo) :: Info(1)
INTEGER N,error
REAL (KIND=qPrec) :: dx,dy,dz
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:) :: ut,vt,wt
! first thing allocate structures for sparse matrixes
INTEGER i,j,k,nx,ny,nz
!INTEGER, PARAMETER :: solver=PARDISO_SOLVER
INTEGER, PARAMETER :: solver=IT_FOM
!INTEGER, PARAMETER :: solver=DCG_SOLVER
LOGICAL iterate








! populate Info structures
Info(1)%mX=(/60,40,14/)

nx=Info(1)%mX(1);ny=Info(1)%mX(2);nz=Info(1)%mX(3)

Info(1)%nDim=3 ! 3D problem
! now define geometry of the problem
Info(1)%xLower=(/10000./2200.,-15000./2200.,0./)

Info(1)%xUpper=(/35000./2200,15000./2200,1./)

Info(1)%dX=(Info(1)%xUpper(1:3)-Info(1)%xLower(1:3))/Info(1)%mX
Info(1)%dmax=0.

! assume all Neumann conditions
Info(1)%mthbc(1:6)=0


!now copy info to elli structure


ALLOCATE(Info(1)%ElliInfo) ! allocate space for ellikit 
Info(1)%ElliInfo%ndim=Info(1)%ndim ! dimension of problem
Info(1)%ElliInfo%mX=Info(1)%mX ! size of grid
Info(1)%ElliInfo%dX=Info(1)%dX ! grid increments
Info(1)%ElliInfo%mthbc=Info(1)%mthbc ! definition of boundary conditions




Call setprob(Info(1))
PRINT*,'Dimension of problem',Info(1)%mX
ALLOCATE(ut(1:ny,1:nz),vt(1:nx,1:nz),wt(1:nx,1:ny))
CALL DEFINE_TOPO(Info(1)) ! define the quantities used to calculate the metric
error=CREATE_OPERATORS(Info(1)%ElliInfo)

! Initially zero in the boundary conditions to be safe
DO i=1,2*Info(1)%nDim
   Info(1)%ElliInfo%bc(i)%p=zero
END DO
! check consistency of BCs
IF(.NOT.check(Info(1)%ElliInfo,CHECK_BC)) THEN
   PRINT*,'BC are not consistent'
   STOP
END IF
Info(1)%ElliInfo%qe=zero
! assign local pointers to edge velocities and fill in with random values
u=>Info(1)%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info(1)%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);
w=>Info(1)%ElliInfo%qe(1:nx,1:ny,1:nz+1,3)




CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
CALL RANDOM_NUMBER(w);
! random values for edge velocities boundary conditions
CALL RANDOM_NUMBER(ut)
CALL RANDOM_NUMBER(vt)
CALL RANDOM_NUMBER(wt)

        

! only need OP= o use IT wth PARDISO precon
i=POISSON_SOLVERS(Info(1)%ElliInfo,INIT,solver,OP=CSR_PRECOND)
ITERATE=.true.;i=1
DO WHILE(ITERATE)
PRINT*,'iteration ',i
i=i+1


Info(1)%ElliInfo%bc(1)%p(1:ny,1:nz)=u(1,:,:)+ut ! u*-utarget
Info(1)%ElliInfo%bc(2)%p(1:ny,1:nz)=u(nx+1,:,:)+ut
Info(1)%ElliInfo%bc(3)%p(1:nx,1:nz)=v(:,1,:)+vt
Info(1)%ElliInfo%bc(4)%p(1:nx,1:nz)=v(:,ny+1,:)+vt

Info(1)%ElliInfo%bc(5)%p(1:nx,1:ny)=w(1:nx,1:ny,1)+wt
Info(1)%ElliInfo%bc(6)%p(1:nx,1:ny)=w(1:nx,1:ny,nz+1)+wt
!!$


PRINT*,'calculates divergence'
   

CALL take_div(Info(1)%ElliInfo,.true.)


div=>Info(1)%ElliInfo%div
PRINT*,'divergence stat'

!!$


PRINT*,'max divergence before = ',maxval(abs(div))

PRINT*,'rms divergence before = ',sqrt(SUM(div**2)/(nx*ny*nz))
IF(solver==DCG_SOLVER) j=POISSON_SOLVERS(Info(1)%ElliInfo,INIT,solver) ! need to reinit solver
j=POISSON_SOLVERS(Info(1)%ElliInfo,SOLVE,solver,OP=CSR_PRECOND)

! subtract gradient
CALL project(Info(1)%ElliInfo,.true.)

PRINT*,'calculates divergence'
10 CALL take_div(Info(1)%ElliInfo,.false.)
div=>Info(1)%ElliInfo%div

PRINT*,'max divergence after = ',maxval(abs(div))

PRINT*,'rms divergence after = ',sqrt(SUM(div**2)/(nx*ny*nz))
ITERATE=sqrt(SUM(div**2)/(nx*ny*nz))>1e-12
ENDDO

i=POISSON_SOLVERS(Info(1)%ElliInfo,FINISH,solver,OP=CSR_PRECOND)
PRINT*,'parison_mopup'
error=FREEUP_OPERATORS(Info(1)%ElliInfo)

DEALLOCATE(Info(1)%topo%h,Info(1)%topo%hx,Info(1)%topo%hy,Info(1)%topo%hxx,Info(1)%topo%hyy)
DEALLOCATE(Info(1)%topo)
DEALLOCATE(ut,vt,wt)
DEALLOCATE(Info(1)%ElliInfo)
STOP

END

